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Abstract 

The rotation rate of the solar radiative zone is an important diagnostic for angular-momentum 
transport in the tachocline and below. In this paper we study the contribution of viscous and mag- 
netic stresses to the global angular-momentum balance. By considering a simple linearized toy 
model, we discuss the effects of field geometry and applied boundary conditions on the predicted 
rotation profile and rotation rate of the radiative interior. We compare these analytical predictions 
with fully nonlinear simulations of the dynamics of the radiative interior, as well as with observa- 
tions. We discuss the implications of these results as constraints on models of the solar interior. 



1. Introduction 

Helioseismic inference of the rotation profile of the solar interior has revealed two spatially- 
distinct regions: an outer differentially -rotating shell surrounding an inner uniformly rotating core 
(Christensen-Dalsgaard & Schou, 1988; Kosovichev, 1988; Brown et al. 1989; Dziembowski et 
al. 1989). The transition between the two regions, the solar tachocline, is located precisely at the 
base of the solar convection zone. It is surprisingly sharp with an average width no larger than 
a few percent of the solar radius (Charbonneau et al. 1999, Elliott & Gough, 1999). Our under- 
standing of this peculiar rotation profile has steadily marched on in the past three decades, bene- 
fiting greatly from the high-performance computing revolution. Today, the field is ripe for more 
quantitative comparisons between models and observations, and has begun focusing on specific 
reference points, such as the overall pole-to-equator difference in the rotation rate, the inclination 
of the isorotation contours, the thickness of the tachocline and finally, the subject of this paper, the 
rotation rate of the radiative interior. 

If one assumes that the solar interior is in a dynamically quasi-steady state, then the rota- 
tion rate of the radiative zone can be thought of as a weighted average of the rotation profile 
observed near the base of the convection zone: 

Qcz(^) — f^eq( 1 - 6(2 COS^ 5 - (34 cos"* 6*) , ( 1 ) 
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where Q.^/2t^ = 463nHz, a2 = 0.17 and = 0.08 (Schou et al. 1998; Gough 2007). Hence we can 
formally write 

/ W (e)flUO) sinOde , (2) 
Jq 

where the weight function W{9) depends uniquely on the nature of angular-momentum transport 
in the tachocline. Observations provide us with relatively precise measurements of firz (Schou et 
al. 1998), with 

Q^Jln = 430nHz ^ O^z 0.93fieq • (3) 
Can this information be used to constrain theoretical models of the solar interior? 

If angular-momentum transport in the radiative zone and the tachocline were purely viscous, 
the rotation rate of the deeper regions (r 0) would be the same as that of the mean specific 
angular momentum of the convection zone (Oilman, Morrow & DeLuca, 1989), or in other words 
WyisciO) = sin^ 6 so that 

a2 3a4 ' 
y~ 35" 

Of course, purely viscous transport cannot account for the observed uniform rotation. Spiegel & 
Zahn (1992) proposed the first tachocline model, in which angular momentum is primarily trans- 
ported by anisotropic turbulence. Since the stratification of the radiative zone strongly inhibits 
radial fluid motions, they modelled the effects of this turbulence as a dominantly horizontal dif- 
fusion process, and found that the interior would indeed relax to uniform rotation beneath the 
tachocline with an angular velocity 



^^visc = ( 1 - ^ - ^ ) Oeq - 0.959Oeq • (4) 



3a2 5fl4 



^^SZ= ( 1-^-^ )aq-0.908neq. (5) 



Gough & Mclntyre (1998, GM98 hereafter) later argued against this model on the ground that 
two-dimensional turbulence does not act to diffuse angular momentum horizontally (e.g. Tobias, 
Diamond & Hughes 2007). Moreover, the predicted value Qsz is sufficiently far from the observed 
value to be ruled out by the observations. 

In the past decade, magnetized models have become more widely accepted as the simplest 
explanation for the observed uniform rotation of the radiative zone (Riidiger & Kitchatinov 1997; 
GM98; see Garaud, 2007 for a review). A large-scale primordial magnetic field, strictly confined 
beneath the convection zone can indeed robustly maintain a state of uniform rotation through Fer- 
raro's law of isorotation (Ferraro, 1937). The confinement of the primordial field is thought to 
result from interactions with large-scale meridional flows originating from the convection zone, as 
originally proposed by GM98, and first shown numerically by Garaud & Garaud (2008) (GG08 
hereafter). 
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However, little attention has been given so far to the predicted rotation rate of the radia- 
tive zone in these models. The linear simulations of Rudiger & Kitchatinov (1997), which as- 
sume a given confined poloidal field structure, suggest that ~ 0.97^7eq5 a value which is much 
larger than observations. The nonlinear simulations of GG08 on the other hand suggest that 
fij-z — 0.87^7eq, a value which is far too low. Can we understand these numerically determined 
values in terms of simple force-balance arguments? We show in this paper that it is indeed pos- 
sible. Moreover, much can be learned from this exercise in terms of relating models to the real 
Sun. 

We begin in §2 by describing a linearized toy model of the radiative zone in which the poloidal 
component of the field is fixed (Rudiger & Kitchatinov 1997; MacGregor & Charbonneau 1999). 
As in the work of MacGregor & Charbonneau (1999), we study two different field geometries: 
an open dipole (§|3]) and a confined dipole (iH). In both cases we study the effect of boundary 
conditions on the predicted rotation profile and rotation rate of the radiative interior. We discuss 
the implications of our findings by comparing these toy-model predictions with fully nonlinear 
simulations of the dynamics of the radiative interior in ^ and conclude in §6. 



2. A toy model of angular momentum transport in the solar radiative zone 

Let us consider the simplest possible setup in which to study the interaction between large- 
scale fields and flows in a spherical shell: the homogeneous magnetized spherical Couette flow. 
We thus consider a spherical shell containing a homogeneous incompressible conducting fluid; the 
outer radius of the shell is while the inner radius is ri. The uniform density of the fluid is p, 
its viscosity v and magnetic diffusivity r]. We model the medium outside the spherical shell as a 
solid, and allow this solid to have various conducting properties. The inner core, for radii r < r-^, is 
assumed to be permeated by currents which maintain a poloidal dipolar magnetic field. 

This setup has been extensively studied in the geophysical literature as a model of the Earth's 
interior and is also used as a basis for studying spherical Couette flow dynamo experiments (e.g. 
Dormy, Cardin & Jauk 1998; Hollerbach, 2000; Dormy, Jauk & Soward 2002; HoUerbach, Canet 
& Foumier, 2007). k has also been used in the solar context by Rudiger & Kitchatinov (1997) 
and MacGregor & Charbonneau (1999). As in these latter papers, we neglect the meridional flows 
entirely. The validity of this approximation is discussed in ^ The flow considered is therefore 
defined in the spherical coordinate system (r, 6*, 0) as 

u = (0,0,rsin^Q(r,^)) . (6) 

As a result of this assumption, the fluid within the spherical shell does not influence the imposed 
poloidal field. 
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We introduce the flux function A(r, 6) such that 

Bp = V X ^^e^ 
\ r sin 6' 



(7) 



The toroidal component of the magnetic field 5^ is generated by the Q-effect induced by the az- 
imuthal flow. For simplicity of notation, we introduce the new variable 5(r, 6) such that 



5 = rsin^5<^ , 
so that the total magnetic field can be written as 



rsin^ \ r 86 dr 



(8) 



(9) 



The dynamics of this reduced system are entirely described by the azimuthal component of 
the momentum equation as well as the azimuthal component of the induction equation. These are 
expressed in the spherical coordinate system as (cf. Rudiger & Kitchatinov 1997): 
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(10) 
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Note that if both diffusion terms are neglected, then the equations reduce to 

Vfi X VA oc Bp ■ VQ = , 
V5 X VA oc Bp ■ V5 = . 



(11) 



(12) 
(13) 



The respective solutions of these equations are simple: Vl and S must be constant on poloidal 
magnetic field lines (which are lines of constant A). Equation (fT2l) is an expression of Ferraro's 
isorotation theorem, while (fT3l) expresses the fact that in this steady-state axisymmetric system the 
azimuthal component of the Lorentz force must be zero. 

The boundary conditions are selected as follows. We consider that the inner boundary is 
rotating uniformly with 

fi(ri,^) = ^]i„, (14) 



while the outer boundary is rotating differentially as 

n{r,,9) = n,M 



(15) 
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The angular velocity is selected, as in Garaud (2002) and GG08, in such a way as to guarantee 
that the total torque applied to the system is zero. Hence, we require that 

dri B B \ 

pur^ sin^ 9— + rsme^-^]smede = . (16) 
or 4% J 

Note that in this steady-state calculation, equation (fT6l) only needs to be applied at one particular 
radius r to be valid everywhere. We now calculate i^in for a variety of poloidal field configurations 
and boundary conditions on the azimuthal magnetic field. 



3. Solution for an open field configuration 

In order to represent an open field configuration, we select 

Air,e>=?lrf-^ . (17) 
2 r 

The normalizing constant is chosen so that this flux function represents the poloidal magnetic field 
Bp which is the exact axisymmetric solution of the equation V^Bp = in the whole space, matches 
to a point dipole at r = and decays as r ^ oo, and finally, which has amplitude Bq on the polar 
axis at radius r = rj. 

Dormy, Cardin & Jault (1998) and Dormy, Jault & Soward (2002) presented the first analytical 
studies of the linear dynamics of magnetized spherical Couette flows. The following analysis is 
analogous to their approach in the limit where meridional flows are neglected, but investigates the 
case of a differentially rotating outer boundary. Following their results we expect the presence of 
two boundary layers near the inner and outer boundaries respectively, as well as an internal shear 
layer along the "last connected field line", as shown in Figure [T] This particular field line (C) 
separates the equatorial region {£) from the polar region (P). 

In the equatorial region {£), every field line (JF) originating from the inner core in the Northern 
hemisphere re-enters the core at a symmetric latitude in the Southern hemisphere. In the limit 
of negligible diffusion, the angular velocity must be constant along (JF) implying that the entire 
equatorial region must rotate with angular velocity fiin. The function S must also be constant along 
(JF), but in addition is antisymmetric with respect to the equator. The only possibility is therefore 
5 = everywhere in {£). 

In the polar region (V), field lines originating from a co-latitude 9i on the inner core extend 
out to co-latitude on the outer boundary (see Figure 1), where 

sin^^i sin^^o 

• (to) 
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Fig. 1. — Unconfined field geometry, showing tlie polar region (V) and the equatorial region (S) 
separated by the last connected field line (C). 

The field line C emerges from the inner core at co-latitude 9c, with 

/ \ 1/2 

sin^c=f^J . (19) 

We now consider the solutions of the problem under various types of boundary conditions for 
the magnetic field. 



3.1. Vanishing toroidal field on the boundaries 

We first consider the case where S = rsinOB^ is required to be zero both at the inner and the 
outer boundary: 

Sirue) = Sir,,e) = (20) 
These can be thought of as "insulating" boundary conditions. 



3.1.1. Analytical solutions 
To study the boundary layer near we introduce the scaled variable 



where the typical boundary layer thickness 6o and its form function f(9) both remain to be deter- 
mined. We assume (and later verify) that for low diffusivities 6o is very small compared with the 
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global scales in the system (ri,ro). Substituting this new variable in the governing equations (flOl) 
and ([TT]) . using the chain rule 

l_d_ ,d_^d__md_ 

and keeping only the lowest order terms in 60 yields 

d's ri6j{e) . , on 

^ = ^0 sm6'cos6'— , 



d^n _ Bo rf6J(6) cos 9 dS 
Combining the two equations yields 



(23) 



d^S dS ^ d^n dn 

provided we define 

*2 = ^p4a„d/(«) = ^. (25) 

5q rf cos 6' 

The solutions to this set of equations which remain bounded as ^ — ^ +00 are 



n,(te) = u'rie) + uji'\9)e-^ , (26) 

where the index "o" denotes that this solution is only valid in the outer boundary layer. The 
integrating functions s\^\9) and io^o\9) are related to one another by the equation 

,a)(^) = _(^l!Lpy^',2,in2 0^a)(^). (27) 
A very similar calculation can be done near ri introducing the scaled variable 



with 

(29) 

yielding similar governing equations (see (l24l) ) and therefore the solutions 



S,iC,9) = sf{9)+s\'\9)e-'^ , 

n,i(:,9) = ujf\9)+ul'\9)e-'^ , (30) 
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where the index "i" denotes that this set of solutions is only valid in the inner boundary layer. The 
integrating functions are related to one another by 

,f\e) = (^^'^ sin' eujl'\e) . od 

Li the case of the boundary conditions selected here, the solution in the bulk of the fluid is 
well-approximated by neglecting any effect of dissipation. This result was formally shown by 
Dormy, Cardin & Jault (1998). There, S and VL are constant along magnetic field lines as discussed 
earlier, so we can write the bulk solution as 

5b(r, 9) = (A(r, 9)) and n^{r, 9) = (A(r, 9)) . (32) 



As ^ and C respectively tend to +oo, the boundary solutions must approach the bulk solution 
smoothly. Hence, 

sf\9) = 5b (A(ri, 9)) and s'-^\9) = 5b (A(ro, 9)) , 
J^\9) = ^b (A(ri, 9)) and lu^^\9) = uj^ {A{r,, 9)) . (33) 

Since the co-latitudes 9i and 9o satisfy by definition A(ri, 9^) =A(ro, 6^0), the above relationships can 
be summarized as 

sf\9d = 5|,°'(^o) and 4^9;) = J^\9,) . (34) 



Finally, we apply the boundary conditions to determine the unknown integration functions 
uniquely. Requiring (fT4l) and (fT5l) implies that 

a;f(^) + u;P\0) = fii„, 
uj^^\9) + uj^J\9)^n,M. (35) 

Requiring 5 = on both boundaries implies 

sf\9)+s'^\9) = 0, 

/^\9)+s^^\9) = 0. (36) 



The set of eight equations contained in (1271 ). (|3T1) . (|34l) . (|34l ) and (|36l) can be solved uniquely 
for the eight integration functions. Of particular interest for the following calculation is the expres- 
sion for (x;^': 

iol'\Oi) = -^^{^cM-^J , (37) 
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where 6i and 6*0 are related by equation (fTSl) . 

We now seek to express as a weighted integral over r2cz(^)- Applying (fT6l) on the inner 
boundary, and using (|28l) and (|37l) yields 

'^^^ sin^ ^^d^ = - ^^u^\e,)^e, = , (38) 



dr Jo 5jm 

since, in the equatorial region on the inner core (^^i > 6c), dVt/dr = 0. Changing variables from 9i 
to 9o, and using the actual expression for f(9) transforms this equation to 

tt/2 

COS sin' 9, [QM - ^in] d^o = , (39) 



f^n. = aq(l-^-^) . (40) 



which implies that 
Using the helioseismically determined values for a2 and a4 yields 

nin = 0.93fieq. (41) 

In the asymptotic limit where this expression is valid (e.g. 5i/ri<^ 1) note that fiin is independent 
of the aspect ratio r^/ of the setup, and can be shown to hold for non-homogeneous fluids as well. 
In fact, it only depends on the imposed differential rotation. 



3.1.2. Numerical solutions 

We obtained numerical solutions to the set of linearized equations (flOl) and (fTT)) with bound- 
ary conditions (fT4l) . (fT5l) and (|20l ) using a numerical algorithm adapted from the one developed by 
Garaud (2001). This algorithm involves the expansion of the governing equations upon a truncated 
set of Chebishev functions in = cos 9, and seeks the solution of the remaining set of ODEs by 
Newton-Raphson relaxation. The linear nature of the equations (in the variables S and Q consid- 
ered) guarantees the immediate convergence of the solutions. 

The analytical results presented above are asymptotically independent of the individual values 
selected for Bq, rj and pu, and depend instead only on their combination through the Hartman num- 
ber H, which is the ratio of the geometric mean of the two diffusion times (viscous and magnetic) 
to the Alfven time. From here on, we define the Hartman number of the simulation as 



(42) 




Fig. 2. — Unconfined field solution for H = 4700 and for S = on both boundaries. In our attempt 
to relate these simulations to the solar radiative zone, the outer solid circle marks for reference the 
outer edge of the Sun (r = Rq), and the dotted line marks the radiative-convective interface. Here, 
ro = Q.IRq and n = 0.35i?Q. 

Figure [2] shows the numerical solution for a fairly large value of the Hartman number, H = 
4700. As expected, the equatorial region is indeed rotating with the same angular velocity as the 
core, and S is zero there. In the polar regions, both Q and S are constant along poloidal magnetic 
field lines, except near the inner and outer boundaries, and in the vicinity of the field line (C). The 
rotation rate of the inner core is found to be ilin = 0.937Qeq in this simulation. This numerical 
solution confirms our analytical results. Moreover, as H ^ +oo, we find that 0.98Qeq (see 

Figure |7]). 



3.2. Conducting boundary conditions 

We now turn to the case of "conducting boundary conditions" for the toroidal field, which 
are the same conditions as those used in the numerical simulations of the radiative zone dynamics 
performed by GG08. The material outside of the fluid shell is assumed to be infinite in extent and 
to have the same diffusivity as the fluid inside the shell. The toroidal field satisfies the equation 

V%-^%- = 0, (43) 
H sm 

for r > To and r < r^, with 5,^ — > as r — and r +oo. The solution to this equation is smoothly 
matched onto the solution within the spherical shell at r = rj and r = (see GG08, for detail). 



As first shown by HoUerbach (2000) the characteristics of the solution are now different from 
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the case of insulating boundary conditions, an effect which is studied in detail by Soward & Dormy 
(2009). Remarkably, the bulk of the fluid is not in Ferraro isorotation. To understand the difference 
on a qualitative basis, note that the magnetic field is merely diffusing out of the boundaries and that 
the continuous generation of toroidal field by the fluid motions within the shell is only compensated 
by this dissipation. In other words, nothing but dissipation limits the growth of the amplitude of the 
field. As a result, the amplitude of S is proportional to l/i], so that the approximation Bp ■ Vfi = 
in the bulk of the fluid is no longer valid (the diffusion term is 0(1) compared with the advection 
term). Note that Bp • V5 = still holds. 

Since the bulk equations are notably more difficult to solve analytically in this case (see 
Soward & Dormy 2009 for detail), we only provide a sample numerical solution in Figure |3l The 
parameters for this simulation are exactly the same as for the case shown in the previous section 
- only the magnetic field boundary conditions differ. The two characteristic features mentioned 
above are clearly illustrated in Figure [3l (1) the amplitude of the S field is orders of magnitude 
larger than before, although S is still constant on poloidal magnetic field lines; (2) the angular ve- 
locity profile deviates from Ferraro isorotation. A strong sub-rotating layer appears just interior to 
the field line (C), and the angular velocity of the inner core = 0.863fieq is much slower than in 
the previous calculation (§3.1). Finally, note that the core velocity found in this case does depend 
on the aspect ratio of the system. The variation of Q-^ with H is shown in Figure |7J 




Fig. 3. — Same as Figure |2lbut for conducting boundary conditions (using (V^B)^ = for r > 
and r < ri). Note how the angular velocity profile is not in a state of Ferraro rotation; also note the 
difference in the amplitude of S. 
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4. Solution for a confined field 

In order to represent a confined magnetic field, we select the flux function A(r, 6) as in the 
work of Riidiger & Kitchatinov (1997) to be 

A(r,^) = 5oy (l--^ysin2^, (44) 

where q is the confinement parameter, which we assume is greater or equal to one. Note that when 
q = I, Br = at To while Bg is non zero. On the other hand, for q > I, Bg also vanishes at r = Tq. 
Also note that the constant Bo now defines the amplitude of the magnetic field at r = (i.e. not on 
the polar axis at r = rj). 

We begin by considering the case of insulating boundary conditions, where 5^ = at and 
Tq. Figure m presents a numerical solution for large Hartman number, for = 1. It clearly shows 
that the bulk of the fluid is rotating uniformly with the same angular velocity as the inner core, and 
that the azimuthal magnetic field is zero in the same region. This bulk solution is expected on the 
same symmetry grounds as the ones invoked in the equatorial region of the unconfined field case. 
It matches smoothly onto the applied boundary conditions at r = through a boundary layer which 
is very thin everywhere except near the polar axis. 




Fig. 4. — Numerical solution in the confined field case for q = I and H = 1340, using insulating 
boundary conditions. Here, ri = 0. 1 and ro = 0.7. The angular velocity is constant in the bulk of the 
"radiative zone", with ^7in = 0.972fieq- 

Following the method used in the case of the unconfined field configuration, we begin by 
rescaling the radial variable near with 

'''' 
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where the typical layer width 5o and the form function f(9) remain to be determined. As before, 
we substitute this new variable for r in equations (flOl ) and (fTTI) and keep the lowest order terms in 
5o only, so that 



+ — ^^t^ + cos6'Kt^ 



2 09 



2 m 



+ : — ^V^ + cos( 



2 89 



If we select 



m 



2 m 
1 



^0 sin2/^6' ' 

then the variables fully separate, and the boundary layer equations become 



(46) 



(47) 



^ dev de) dfi[ 



^ de V de 



5/i 



(48) 



with the introduction of the variable /i = cos 9. 



We validate our boundary layer approximation by comparing, in Figure[5l the predicted latitu- 
dinal structure of the boundary layer for (^r = 1 with numerical solutions. The agreement is excellent, 
and similar comparisons for other values of q are also in excellent agreement with (1471 ). 

Even when q=l, solving (l48l) analytically is not entirely triviail]. Luckily, estimating the value 
of the interior angular velocity does not require knowledge of the full boundary layer solution. We 
can obtain a first-order accurate approximation of the derivative dVl / dr across the boundary layer 
with 

ficz(^)-^in 



dn 

dr 



(49) 



for large Hartman number. Applying (fT6l) at and using the approximation (|49l) . we find that the 
angular velocity of the interior is determined by 



which implies 



a, 



/ (l-//2)i+i/'?[a,(/i)-Q,„]d/i = 



5 + 2/^ (7 + 2/^)(5 + 2/^) 



(50) 



(51) 



' An approximate solution can be found for q=\. 
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atitude 

Fig. 5. — Close up of the angular velocity profile near the outer boundary layer for ^ = 1 and 
H - 1340 (see FigureH]). Superimposed on this numerical solution we show the predicted boundary 
layer thickness 5of{6) divided by three for the same parameter values - the two are in excellent 
agreement. 

Note that as ^ ^ oo (when the field is more and more confined), the solution recovers the purely 
viscous case as expected. For the helioseismically determined values of and 04 we get 

l^in ~0.972Qeqfor^ = 1 , 

17in ~ 0.966Qeq for ^ = 2 , (52) 

and so forth. To verify our analysis, we compare the asymptotic values of Vtu^ calculated in equation 
(ISTI) with the numerical solutions for different values of q in Fig. |6l The agreement, for high 
values of the Hartman number, is again excellent. Finally, it can be shown that by contrast with 
the unconfined field case, changing boundary conditions for the toroidal field does not yield a 
different asymptotic answer for VL^^. This can be attributed to the fact that the bulk solution, and 
the geometry of the boundary layer, are the same regardless of the boundary conditions. 



5. Discussion 

In the two previous sections, we studied the axially symmetric, linearized dynamics of ho- 
mogeneous magnetized spherical Couette flows for various geometries of the imposed poloidal 
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Fig. 6. — Calculated core velocity as a function of Hartman number. This figure compares numeri- 
cal simulations and analytical predictions for l^in/ figq for two values of the confinement parameter 
q, for insulating boundary conditions. The diamonds show numerical results for ^ = 1, and the plus 
symbols for q = 2. The dotted line marks the analytical asymptotic limit for ^ = 1 and the dashed 
line for ^ = 2 (see equation (eq:ominconf2)). In the low Hartman number limit on the other hand, 
all curves converge to the value f^visc/^eq (solid line). 



field and for various types of magnetic boundary conditions. We now discuss how these simplified 
models may help us make sense of the dynamics of the solar radiative zone and the tachocline. 



5.1. Comparison with the numerical simulations of GG08 

Recently, GG08 presented a set of numerical simulations of the solar radiative zone and the 
tachocline based on the Gough & Mclntyre model (GM98). By contrast with the toy model studied 
above, GG08 attempt to model the dynamics of the radiative interior as accurately as possible and 



- 16- 



solve the following system of equations 



pu- Vu+2pQo X u = -Vp-pV$+j X B+/j,V- n , 
V-(pC) = 0, 

pru-V5 = v-(Awf), 



p p T 



V X (u X B) = V X (f^f]V X B) , 
VB = 0, 



(53) 



for the velocity field in the rotating frame u = (Ur, ug, r sin 6Vl), for the magnetic field B = {B, ,B0,B^) 
and for the density (p), temperature {T) and pressure (p) perturbations. In these equations, VLq is 
equal to l^visc (see equation dH)), the background stratification for p (the density), f (the temper- 
ature), s (the entropy), $ (the gravitational potential) and p (the pressure) are given by Model S 
of Christensen-Dalsgaard et al. (1996), and the diffusion coefficients v (for the viscosity), fj (for 
the magnetic diffusivity) and k (for the thermal conductivity) are calculated according to Gough 
(2007) (see also GG08). Large factors /j,, and fk multiply these respective diffusivities to help 
numerical convergence. The boundary conditions used on the inner core are: 

• impermeable and no-slip, with Q = - VIq and where is deduced from (fT6l) . 

• electrically conducting for the electric/magnetic field (e.g. the magnetic field satisfies V^B = 
in the inner core, matches on to a point dipole at r = 0, and matches continuously to the 
solution in the spherical shell at r = ri), 

• thermally conducting for the temperature field (e.g. T satisfies V^T = in the inner core, 
and smoothly matches on to the solution in the shell at r = n). 

The boundary conditions used on the outer boundary are: 

• the velocity field matches smoothly onto an imposed velocity field 



• electrically conducting for the electric/magnetic field (e.g. the magnetic field satisfies V^B = 
for r > To, vanishes as r — > oo and smoothly matches onto the shell solution at the boundary, 

• thermally conducting for the temperature field (e.g. T satisfies V^T = for r > r^, and 
smoothly matches on to the solution in the shell at r^. 




M'^Me^roSin^(Qcz-^o)) , 
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The numerical solutions of these equations and boundary conditions, as computed by GG08, 
exhibit many of the dynamical properties of the tachocline first discussed by GM98. In order to 
prevent the propagation of the convection zone shear into the interior (as illustrated in Figures 
[21 and |3]), confining the primordial field to the radiative zone appears to be necessary. GM98 
argued that large-scale meridional flows downwelling from the convection zone would naturally 
interact nonlinearly with the underlying field and confine it beneath the bulk of the tachocline. The 
convection zone flows are by the same mechanism prevented from penetrating more deeply into 
the radiative zone, thereby satisfying a variety of observational constraints on chemical mixing 
near the base of the convection zone (e.g. Elliott & Gough 1999). The tachocline would thus be 
a well-ventilated region, spatially separated from the magnetically-dominated interior by a very 
thin advection-diffusion layer. All of these features are qualitatively well- accounted for in the 
simulations of GG08. 

However, the angular velocity of the bulk of the radiative interior in the simulations of GG08 
is much lower than the observed value: Qin — 0.87f2eq in the limit of large Hartman number. 
Surprisingly, the same value is found for a fairly wide range of assumed convection zone flow 
amplitudes and profiles uf'{6). Given that simulations robustly insist on selecting this particular 
interior angular velocity, can we understand it in terms of the simple dynamics studied in the toy 
model? Surprisingly, it appears that we can. 

To see this quantitatively. Figure |7] compares the core angular velocity predictions for the 
stratified, nonlinear calculations of GG08 with our toy-model calculations for an open dipole in a 
similar experimental setup (i.e. the same aspect ratio and with conducting boundary conditions). 
The agreement between the two sets of simulations as a function of Hartman number is quite 
remarkable, in spite of the over-simplified nature of the toy model. 

We did in fact expect the agreement to be very good for high values of the diffusivities (low 
values of H). Indeed, since the background is strongly stratified, only very slow meridional flows 
can penetrate into the radiative interior (see Garaud & Brummell 2008). The magnetic Reynolds 
number of these flows is therefore also low, and they fail to have any influence on the poloidal 
field. The magnetic field then naturally relaxes to its fundamental eigenmode, which is the open 
dipole. These combined factors together imply that the linearized problem studied in §|3] should 
be (and is indeed found to be) a good approximation to the overall angular momentum balance of 
the system. It correctly predicts the asymptotic limit // ^ 0, as well as the somewhat surprising 
bifurcation around H =1. It is also interesting to note that radial variations in p, v and ri do not 
appear to affect the predictions for the interior angular velocity much. This can be formally shown 
in the case of insulating boundary conditions but not in the case of conducting boundary condi- 
tions. Nevertheless, it appears that the agreement approximately holds for conducting boundary 
conditions as well. 
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Fig. 7. — Calculated angular velocity of the core as a function of Hartman number. This fig- 
ure compares the linear solutions in the toy-model (lines) and fully nonlinear numerical solutions 
(symbols) for two different types of boundary conditions. The conducting boundary conditions, 
shown in the dotted line and diamond symbols, were also used by GG08. In the limit of large 
Hartman number, the interior rotation rate i^in tends to about O.SVQeq- In the second case, shown 
with the dashed line and triangles, we changed the lower boundary to be insulating. The asymp- 
totic limit in both toy model and in the full numerical simulations then appears to be closer to the 
observations, with Clin 0.93fieq- 

The good agreement between the predicted core angular velocities in the toy model and in 
GGOS's simulations, for low values of the diffusivities (high values of H), is much more surpris- 
ing. The lowest-diffusivities simulations presented by GG08, which correspond to the right-most 
symbols in Figure |7] (see also Figure (8]), have a magnetic field geometry which deviates signifi- 
cantly from the purely dipolar open-field configuration studied in O Moreover, the same simula- 
tions clearly show the existence of a region where angular-momentum transport is operated by the 
meridional flows rather than by the magnetic field, as predicted by GM98. Both phenomena are a 
natural consequence of the increasingly nonlinear nature of the interaction between the primordial 
field and the assumed convection zone flows as the diffusivities are lowered (equivalently, as H is 
increased); both should by and large invalidate the applicability of the toy model. Nevertheless, 
even then we find that the open-dipole toy model adequately predicts the core angular velocity of 
the fully nonlinear numerical simulations, for the parameter values considered. 
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Whether this statement would continue to hold for even lower values of the diffusivities (even 
higher values of H) is a priori unlikel>@. As the importance of angular-momentum transport by 
the meridional flows increases in relation to viscous transport, we expect significant deviations 
away from the toy model predictions to occur. Indeed, while GG08 failed to achieve low enough 
diffusion parameters in their simulations to test this hypothesis, they also presented another set of 
simulations for artificially high convection zone flow amplitudes (see their Figure 1 1), in which the 
calculated core velocities do deviate significantly away from the toy-model predictions. Moreover, 
one could also expect that as the magnetic field becomes even more confined to the interior, the 
open-dipole configuration will lose relevance in favor of the closed-dipole configuration. Since 
we have shown that the closed-dipole predictions are closer to l^in ~ O.QV^eqj we may expect the 
predicted core velocity in the full simulations to increase towards this value as the diffusivities are 
decreaseci^. 



5.2. Sensitivity to boundary conditions 

The core velocity found in the linearized model is sensitively dependent on the assumed mag- 
netic field boundary conditions. This was shown in ^ where changing the boundary conditions 
from insulating to conducting had a profound impact on the nature of the solution. In fact, it can 
also be shown that the same happens by changing from conducting conditions to having even just 
one boundary condition where B^ = (see below). But far more importantly, this sensitive depen- 
dence on boundary conditions is also found in the numerical solutions of the full set of equations 
(El. 

We ran a separate set of simulations for both the linearized toy model and for the full nonlinear 
model, where the outer boundary is again conducting but the inner core is now insulating (in 
the sense that 5^ = at r = rj). The predictions for the core angular velocity as a function of 
the Hartman number are shown in Figure |7J We found that in the large-Hartman number limit, 
the "asymptotic" value of f^in in the full nonlinear model is much closer to the observed value, 
0.93Qeqj as also predicted by the toy model. Figure [8] compares the numerical results of the full 
simulations with the different sets of boundary conditions with one-another. While the streamlines 
and the poloidal field lines (as well as the temperature and density profiles, not shown) appear to 
be relatively unaffected by the new boundary condition at r^ the angular velocity and toroidal field 
profiles are notably different. Following the trends of the toy model, the amplitude of the toroidal 
field is significantly lower and the rotation profile is much closer to Ferraro isorotation when at 



^although cannot be ruled out. 

^unless angular-momentum transport by the flows acts just in the opposite way, see footnote 2. 



-20- 



least one of the boundary conditions is not conducting. 
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Fig. 8. — Comparison between the numerical solutions of (l53l) with two different sets of boundary 
conditions on the toroidal field. Top: Conducting boundary conditions everywhere. Bottom: The 
top boundary is conducting while the bottom boundary is insulating. All other simulation parame- 
ters are exactly the same for the two runs: Bq = 7T, f^ = Sx 10^, = 8 x 10^ and = 8 x 10^. The 
strip beneath each quadrant zooms into the region near the outer boundary, for r E [O.65,O.7]i?0. 
The numbers represent latitude. In the streamlines panel, solid lines denote clockwise flows and 
dotted lines anti-clockwise flows. 

The implications of these findings are quite important. Short of simulating the entire solar 
interior including the turbulent convection zone and its effect on magnetic fields, one needs to make 
assumptions on the nature of the radiative-convective interface. It appears that models in which 
the fluid shell is contained in a conducting solid of infinite extent are somewhat pathological in 
nature, as they allow an unphysically high magnetic field amplitude to build up thus breaking away 
from Ferraro isorotation (Soward & Dormy 2009). Meanwhile, insulating boundary conditions 
seem to be an a priori equally poor physical representation of both the inner core and of the 
radiative-convective interface. In reality, one may either expect other physical phenomena to limit 
the toroidal field amplitude within the radiative zone (e.g. magnetic instabilities), or at the very 
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least note that the Sun is not infinite in extent, so that 5,^ = in the vacuum outside of Rq. So, as 
strange as it may seem, having at least one insulating boundary may actually be more physically 
realistic than the boundary conditions originally used by GG08. 



5.3. Implications for models and observations of the solar interior 

Aside from the demonstrably odd case of the conducting boundary conditions described in 
the previous section, all simple analytical models of the radiative zone presented so far predict 
an angular velocity fiin close to the observations: in all cases (see Table 1), 0.908fieq < < 
0.972r2eq- Crucially, all of the estimates presented in Table 1 are independent of p, u, r], Bq and 
of the aspect ratio ro/rj, in the asymptotic limit of large Hartman number. In addition, Garaud 
(2002) also studied the hydrodynamic case of spherical Couette flow between one differentially 
and one uniformly rotating sphere, and showed that when the gap width is of the order of the 
observed thickness of the tachocline, the predicted angular velocity of the interior is also, perhaps 
coincidentally, close to the observed value of fir^ = 0.93fieq- 

Table 1 has a somewhat ironic property: the models which are a priori the most unphysical, 
or the poorest representation of the solar interior are the ones which actually seem to fare the best 
in terms of predicting Qjn close to the observed value. Indeed, recall that the open-dipole case has 
a non-uniformly rotating radiative zone, while the hydrodynamic spherical Couette flow (Garaud, 
2002) assumes the fluid to be confined between two impermeable spherical shells. 

There are several lessons to be learned from this work. As mentioned in §1, all predictions 
for necessarily involve a weighted integral over i7cz(^)- Moreover, the spherical geometry of 
the problem implies that the weight function is typically biased towards the equatorial regions 
- in other words, Vlin is more sensitive to a2 than to a4. As a result, we see that the spread in 
predictions for Qjn is relatively small, and one should neither be surprised to see many different 
models predicting similar values, nor that some should lie coincidentally close to the observed one. 

Nevertheless, it is equally interesting to see that fiin in the closed-dipole model, which is per- 
haps the "closest" (in relative terms) to what one may expect from the tachocline dynamics, is sig- 
nificantly different from the observations. This implies one of two things: either meridional flows 
(or perhaps anisotropic turbulent stresses) are a non-negligible contribution to angular-momentum 
transport in the tachocline or (if they are negligible) the true angular velocity profile near the base 
of the convection zone deviates significantly away from the one used here (see equation ([T])). He- 
lioseismology may be able to help distinguish between these two alternatives. 
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6. Conclusion 

We have studied, analytically and numerically, the predicted angular velocity profile of the 
solar radiative zone under various model assumptions. Our overall conclusions have implications 
for future modeling, and implications for future observations. 

In terms of modeling, we have illustrated how crucial the selection of magnetic boundary 
conditions can be to the calculated solution, an effect which has only recently been fully appreci- 
ated (see the detailed study by Soward & Dormy, 2009). Assuming, as previous models have done 
(Garaud, 2002; Brun & Zahn 2006; GG08), that the radiative zone is contained within a homoge- 
neous conducting medium of infinite extent allows unphysically large toroidal field amplitude to 
build up. This case is a somewhat pathological limit, since if the toroidal field is somehow forced 
to be zero at a finite radius (e.g. the solar photosphere), or if other mechanisms act to limit its 
amplitude, then the problem does not arise. Nevertheless, how to best represent the presence of the 
solar convection zone remains to be determined. 

In terms of observations, our various calculations have quantified the sensitivity of the angular 
velocity of the interior to the model assumptions: aside from a few exceptional cases which can be 
ruled out (see above) the predicted angular velocity lies roughly in the interval [0.911)eq,0.97Qeq]- 
Angular-momentum balance between viscous stresses and magnetic stresses for a closed-dipole 
suggests that ~ 0.97Qeq- If helioseismic observations can rule out this value entirely, then 
we can conclude from this study that the tachocline is the seat of additional mixing, either in the 
form of large-scale meridional flows (Gough & Mclntyre, 1998), or in the form of small-scale 
turbulence. Although chemical evidence for additional mixing in the tachocline has already been 
put forward (Gough & Mclntyre, 1998; Elliott & Gough, 1999, Rudiger & Pipin, 2001), our work 
provides the first dynamical evidence to this effect. 
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Table 1: Summary of analytical model predictions 



Model type 




eq 


^rz/^eq 


Viscous model'' 


1 12 
^ 5 


3t/4 


0.959 


Anisotropic viscosity'^ 


1 3^2 

7 


504 
21 


0.908 


Open dipole field'^ 


1 «2 
^ 3 


^4 

6 


0.930 


Confined dipole field'^ 1 - 


02 3(14 
5+2/9 (5+2/9)(7+2/9) 


0.959 - 0.972 



"Using fl2 = 0.17 and 04 = 0.08, Schou et al. (1998) 
''Oilman, Morrow & DeLuca (1989) 
"Spiegel &Zahn (1992) 
''see ED 
"see g] 



